/**
 * @author Gao Jian
 * @email gaoj@njupt.edu.cn
 * @create date 2021-06-03 10:23:40
 * @desc
 * 关于tiff掩膜数据的处理，包括连接破碎的小块成大块，删除块内空洞，较小的块
 */

#define cimg_use_tiff

#include <CImg.h>
#include <gdal/gdal.h>
#include <iostream>

using namespace std;
using namespace cimg_library;
void ain_regions_cpu(unsigned char *data, unsigned char value, unsigned sz);
void bin_regions_cpu(unsigned char *data, unsigned char value, unsigned sz);
void dilateparts_cpu(unsigned char *output, unsigned char *input,
                     unsigned char value, unsigned int nw, unsigned int nh);
int clean_parts_cpu(unsigned char *data, unsigned char val1, unsigned char val2,
                    unsigned int threshold, unsigned int nw, unsigned int nh,
                    const unsigned verbose = 0);

int main(int argc, char **argv) {
  const char *fi = cimg_option("-f", (char *)(NULL), "TiFF File Name");
  const char disp = cimg_option("-dis", 0, "Visualize data");
  const char *pmode = cimg_option("-p", "D9F0E9C0", "Procedure mode, 'DFEC'");
  const unsigned char vala = cimg_option("-va", 13, "binarize data value");
  const unsigned char valb = cimg_option("-vb", 30, "binarize data value");
  const unsigned char vald = cimg_option("-vd", 0x01, "data value");
  const unsigned char valn = cimg_option("-vn", 0x00, "nonvalid value");
  const unsigned char jmode = cimg_option("-jm", 0xFF, "Joing operation mode");
  const unsigned resl = cimg_option("-r", 0, "resolution (meter).");
  const unsigned thf = cimg_option("-tf", 100, "threshold for filling");
  const unsigned thc = cimg_option("-tc", 100, "threshold for cleaning");
  const unsigned verbose = cimg_option("-v", 0, "verbose.");
  const unsigned device_id = cimg_option("-d", 0, "device ID.");

  if (verbose)
    printf("========== Dilate & Fill & Erode & Clean raster parts. "
           "==========\n");

  if (!fi)
    return -1;

  // =========================================================================
  GDALDatasetH hDataset;
  GDALAllRegister();
  hDataset = GDALOpen(fi, GA_Update);
  if (hDataset == NULL) {
    return -1;
  }

  double adfGeoTransform[6];

  GDALGetGeoTransform(hDataset, adfGeoTransform);
  unsigned resolution = resl;
  double arc = 6378245. * 3.14159265358 * 2 / 360;
  if (resl == 0) {
    if (strstr(GDALGetProjectionRef(hDataset), "degree"))
      resolution =
          unsigned((adfGeoTransform[1] - adfGeoTransform[5]) / 2 * arc + 0.5);
    else
      resolution =
          unsigned((adfGeoTransform[1] - adfGeoTransform[5]) / 2 + 0.5);
  }

  GDALRasterBandH hBand;
  // int nw, nh;
  hBand = GDALGetRasterBand(hDataset, 1);
  int nw = GDALGetRasterBandXSize(hBand);
  int nh = GDALGetRasterBandYSize(hBand);
  int sz = nw * nh;
  // GDALGetBlockSize(hBand, &nw, &nh);
  // printf(
  //     "Block=%dx%d Type=%s, ColorInterp=%s\n", nw, nh,
  //     GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
  //     GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
  CImg<unsigned char> datain(nw, nh);
  CPLErr error = GDALRasterIO(hBand, GF_Read, 0, 0, nw, nh, datain.data(), nw,
                              nh, GDT_Byte, 0, 0);

  // =========================================================================
  if (verbose)
    printf("========== Labelling Built Up Area (%dm). ==========\n",
           resolution);

  CImg<unsigned char> odat = datain;

  CImg<unsigned char> tmp = odat;

  CImg<unsigned char> dataout(nw, nh, 1, 1, 0);
  dataout = datain;

  unsigned char nn = strlen(pmode) / 2;

  for (char i = 0; i < nn; i++) {

    const char *proc = pmode + i * 2;
    char op[2];
    op[1] = '\0';

    if (*proc == 'A' || *proc == 'a') { // 二值化
      if (verbose)
        printf("  Step [%d/%d]: BINARIZING [%d] parts.\n", i + 1, nn, vala);
      ain_regions_cpu(dataout.data(), vala, sz);
    }

    if (*proc == 'B' || *proc == 'b') { // 二值化
      if (verbose)
        printf("  Step [%d/%d]: BINARIZING [%d] parts.\n", i + 1, nn, valb);
      bin_regions_cpu(dataout.data(), valb, sz);
    }

    if (*proc == 'D' || *proc == 'd') { // 膨胀碎片
      op[0] = proc[1];
      unsigned dnit = atoi(op);
      if (verbose)
        printf("  Step [%d/%d]: DILATING [%d] parts for %d times.\n", i + 1, nn,
               vald, dnit);
      for (int i = 0; i < dnit; i++) {
        dilateparts_cpu(tmp.data(), dataout.data(), vald, nw, nh);
        dataout = tmp;
      }
    }

    if (*proc == 'E' || *proc == 'e') { // 腐蚀碎片
      op[0] = proc[1];
      unsigned enit = atoi(op);
      if (verbose)
        printf("  Step [%d/%d]: ERODING [%d] parts for %d times.\n", i + 1, nn,
               vald, enit);
      for (int i = 0; i < enit; i++) {
        dilateparts_cpu(tmp.data(), dataout.data(), valn, nw, nh);
        dataout = tmp;
      }
    }

    if (*proc == 'F' || *proc == 'f') { // 填充小洞
      if (verbose)
        printf("  Step [%d/%d]: FILLING [%d] parts with threshold [%d m^2].\n",
               i + 1, nn, valn, thf);
      clean_parts_cpu(dataout.data(), valn, vald,
                      thf / (resolution * resolution), nw, nh);
    }

    if (*proc == 'C' || *proc == 'c') { // 删除小片
      if (verbose)
        printf("  Step [%d/%d]: CLEANING [%d] parts with threshold [%d m^2].\n",
               i + 1, nn, vald, thc);
      clean_parts_cpu(dataout.data(), vald, valn,
                      thc / (resolution * resolution), nw, nh);
    }
  }

  // =========================================================================
  error = GDALRasterIO(hBand, GF_Write, 0, 0, nw, nh, dataout.data(), nw, nh,
                       GDT_Byte, 0, 0);

  GDALClose(hDataset);

  // =========================================================================

  if (verbose)
    printf("========================= <<<< Finished. >>>> "
           "=========================\n");
  // if (fo)
  //   dataout.save(fo);

  if (disp)
    (odat, dataout).display("INPUT | OUTPUT");

  return 0;
}
